
% This script performs the Boostrap Monte Carlo simulations

clc;
clear;

tic

% Get original estimates
regression_output = readtable('RegressionOutput.csv');

regression_output = table2array(regression_output);

Y = -regression_output(:,1:3)';

Sigma = [regression_output(:,4), regression_output(:,7), regression_output(:,8);...
         regression_output(:,7), regression_output(:,5), regression_output(:,9);...
         regression_output(:,8), regression_output(:,9), regression_output(:,6)];

% Get bootstrap estimates
simulation_output = readtable('SimulationOutput.csv');

simulation_output = table2array(simulation_output);

num = height(simulation_output);

TO = -simulation_output(:,1);
CO = -simulation_output(:,2);
B = -simulation_output(:,3);
Var_TO = simulation_output(:,4);
Var_CO = simulation_output(:,5);
Var_B = simulation_output(:,6);
Cov_TO_CO = simulation_output(:,7);
Cov_TO_B = simulation_output(:,8);
Cov_CO_B = simulation_output(:,9);

beta0 = Y;

CO = CO - beta0(2);

beta0(2) = 0;

alpha = 0.05;

bounded = [1;1;0];

bounds = [0;0;0];

R = [1 0 0;0 1 0;-1 -1 1];

cover_SSCI_TO = zeros(num,1);
cover_SSCI_CO = zeros(num,1);
cover_SSCI_B = zeros(num,1);
cover_SSCI_I = zeros(num,1);
cover_SSCI_B0 = zeros(num,1);
cover_SSCI_I0 = zeros(num,1);

cover_SCI_TO = zeros(num,1);
cover_SCI_CO = zeros(num,1);
cover_SCI_B = zeros(num,1);
cover_SCI_I = zeros(num,1);

excess_length_SSCI_TO = zeros(num,1);
excess_length_SSCI_CO = zeros(num,1);
length_SSCI_B = zeros(num,1);
length_SSCI_I = zeros(num,1);
length_SSCI_B0 = zeros(num,1);
length_SSCI_I0 = zeros(num,1);

excess_length_SCI_TO = zeros(num,1);
excess_length_SCI_CO = zeros(num,1);
length_SCI_B = zeros(num,1);
length_SCI_I = zeros(num,1);

ratio_excess_length_TO = zeros(num,1);
ratio_excess_length_CO = zeros(num,1);
ratio_length_B = zeros(num,1);
ratio_length_I = zeros(num,1);
ratio_length_B0 = zeros(num,1);
ratio_length_I0 = zeros(num,1);

for it=1:num
    % TO, CO, and B
    theta_b = [TO(it);CO(it);B(it)];
    
    Sigma_b = [Var_TO(it),    Cov_TO_CO(it), Cov_TO_B(it);...
               Cov_TO_CO(it), Var_CO(it),    Cov_CO_B(it);...
               Cov_TO_B(it),  Cov_CO_B(it),  Var_B(it)];
           
    std_errors_b = sqrt(diag(Sigma_b));
    
    % TO
    SSCI_TO = SSCI(theta_b,Sigma_b,1,1,alpha,bounded,bounds);
    
    if SSCI_TO(1) < beta0(1)
        cover_SSCI_TO(it) = 1;
    end
    
    if theta_b(1)-norminv(1-alpha)*std_errors_b(1) < beta0(1) 
        cover_SCI_TO(it) = 1;
    end 

    excess_length_SSCI_TO(it) = theta_b(1)-SSCI_TO(1);

    excess_length_SCI_TO(it) = norminv(1-alpha)*std_errors_b(1);
    
    ratio_excess_length_TO(it) = (theta_b(1)-SSCI_TO(1))/(norminv(1-alpha)*std_errors_b(1));

    % CO
    SSCI_CO = SSCI(theta_b,Sigma_b,2,1,alpha,bounded,bounds);
    
    if SSCI_CO(1) < beta0(2)
        cover_SSCI_CO(it) = 1;
    end
    
    if theta_b(2)-norminv(1-alpha)*std_errors_b(2) < beta0(2) 
        cover_SCI_CO(it) = 1;
    end 
    
    excess_length_SSCI_CO(it) = theta_b(2)-SSCI_CO(1);

    excess_length_SCI_CO(it) = norminv(1-alpha)*std_errors_b(2);
    
    ratio_excess_length_CO(it) = (theta_b(2)-SSCI_CO(1))/(norminv(1-alpha)*std_errors_b(2));
    
    % B
    SSCI_B = SSCI(theta_b,Sigma_b,3,0,alpha,bounded,bounds);
    
    if ~isempty(SSCI_B) && SSCI_B(1) < beta0(3) && SSCI_B(2) > beta0(3)
        cover_SSCI_B(it) = 1;
    end
    
    SSCI_B0 = SSCI(theta_b-[beta0(1);0;0],Sigma_b,3,0,alpha,bounded,bounds);
    
    if ~isempty(SSCI_B0) && SSCI_B0(1) < beta0(3) && SSCI_B0(2) > beta0(3)
        cover_SSCI_B0(it) = 1;
    end
    
    SCI_B = [theta_b(3)-norminv(1-alpha/2)*std_errors_b(3),theta_b(3)+norminv(1-alpha/2)*std_errors_b(3)];
    
    if SCI_B(1) < beta0(3) && SCI_B(2) > beta0(3)
        cover_SCI_B(it) = 1;
    end
    
    if ~isempty(SSCI_B)
        length_SSCI_B(it) = SSCI_B(2)-SSCI_B(1);
    else
        length_SSCI_B(it) = 0;
    end
    
    length_SCI_B(it) = SCI_B(2)-SCI_B(1);
    
    ratio_length_B(it) = length_SSCI_B(it)/length_SCI_B(it);
    
    if ~isempty(SSCI_B0)
        length_SSCI_B0(it) = SSCI_B0(2)-SSCI_B0(1);
    else
        length_SSCI_B0(it) = 0;
    end
    
    ratio_length_B0(it) = length_SSCI_B0(it)/length_SCI_B(it);
    
    % I
    theta_b = R*theta_b;

    Sigma_b = R*Sigma_b*R';
    
    std_errors_b = sqrt(diag(Sigma_b));
    
    SSCI_I = SSCI(theta_b,Sigma_b,3,0,alpha,bounded,bounds);
    
    if ~isempty(SSCI_I) && SSCI_I(1) < R(3,:)*beta0 && SSCI_I(2) > R(3,:)*beta0
        cover_SSCI_I(it) = 1;
    end
    
    SSCI_I0 = SSCI(theta_b-R*[beta0(1);0;0],Sigma_b,3,0,alpha,bounded,bounds);
    
    if ~isempty(SSCI_I0) && SSCI_I0(1) < (R(3,:)*(beta0-[beta0(1);0;0])) && SSCI_I0(2) > (R(3,:)*(beta0-[beta0(1);0;0]))
        cover_SSCI_I0(it) = 1;
    end
    
    SCI_I = [theta_b(3)-norminv(1-alpha/2)*std_errors_b(3),theta_b(3)+norminv(1-alpha/2)*std_errors_b(3)];
    
    if SCI_I(1) < R(3,:)*beta0 && SCI_I(2) > R(3,:)*beta0
        cover_SCI_I(it) = 1;
    end
    
    if ~isempty(SSCI_I)
        length_SSCI_I(it) = SSCI_I(2)-SSCI_I(1);
    else
        length_SSCI_I(it) = 0;
    end
    
    length_SCI_I(it) = SCI_I(2)-SCI_I(1);
    
    ratio_length_I(it) = length_SSCI_I(it)/length_SCI_I(it);
    
    if ~isempty(SSCI_I0)
        length_SSCI_I0(it) = SSCI_I0(2)-SSCI_I0(1);
    else
        length_SSCI_I0(it) = 0;
    end
    
    ratio_length_I0(it) = length_SSCI_I0(it)/length_SCI_I(it);
    
end

[mean(cover_SSCI_TO)*100 mean(cover_SSCI_CO)*100 mean(cover_SSCI_B)*100 mean(cover_SSCI_I)*100 mean(cover_SSCI_B0)*100 mean(cover_SSCI_I0)*100;...
 mean(cover_SCI_TO)*100 mean(cover_SCI_CO)*100 mean(cover_SCI_B)*100 mean(cover_SCI_I)*100 mean(cover_SCI_B)*100 mean(cover_SCI_I)*100;...
 mean(excess_length_SSCI_TO) mean(excess_length_SSCI_CO) mean(length_SSCI_B) mean(length_SSCI_I) mean(length_SSCI_B0) mean(length_SSCI_I0);...
 mean(excess_length_SCI_TO) mean(excess_length_SCI_CO) mean(length_SCI_B) mean(length_SCI_I) mean(length_SCI_B) mean(length_SCI_I);...
 mean(excess_length_SSCI_TO)/mean(excess_length_SCI_TO) mean(excess_length_SSCI_CO)/mean(excess_length_SCI_CO)...
 mean(length_SSCI_B)/mean(length_SCI_B) mean(length_SSCI_I)/mean(length_SCI_I) mean(length_SSCI_B0)/mean(length_SCI_B) mean(length_SSCI_I0)/mean(length_SCI_I)]

toc

